Maximum Valency Lattice Gas Models 

Srikanth Sastry/ Emilia La Nave,^ and Francesco Sciortino^ 

^ Jawaharlal Nehru Centre for Advanced Scientific Research, 
Jakkur Campus, Bangalore 560064, INDIA 
^ Dipartimento di Fisica and INFM-CRS-SOFT, 
Universitd di Roma La Sapienza, P.le A. Moro 2, 00185 Roma, Italy 

Abstract 

We study lattice gas models with the imposition of a constraint on the maximum number of 
bonds (nearest neighbor interactions) that particles can participate in. The critical parameters, 
as well as the coexistence region are studied using the mean field approximation and the Bethe- 
Peierls approximation. We find that the reduction of the number of interactions suppresses the 
temperature-density region where the liquid and gas phases coexist. We confirm these results 
from simulations using the histogram reweighting method employing grand Canonical Monte Carlo 
simulations. 
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I. INTRODUCTION 



The study of theggl state of matter in colloidal systems is receiving significant attention 
in the last years^, UO, 0]. At the heart of this interest lies the hope that these studies will 
help understanding differences and similarities in the processes of formation of arrested states 
of matter at low packing fraction 0, and in particular the inter-relations between gels and 
glasses. A similarly ambitious additional goal is to provide a new route for understanding the 
intrinsic properties of the process of formation of physical gels in systems more complicated 
than colloids, as the case of gels formed by reversible cross- linking of polymeric chains [sj and 
gelation in protein solutions j^. The full comprehension of both these aggregation processes 
are of extreme importance in the food industry, in the protein crystallization process and in 
the design of novel biomaterialsj^]- 

One of the key questions regards the interplay between the process of formation of a 
long-living network and the process of phase separation in colloid-rich and colloid-poor 
regions^, y, Q, Q. Indeed, the increase of the bond interaction strength (relative to 
the thermal energy) controls both the increase of the inter-particle bond lifetime (and as a 
consequence the lifetime of the spanning network) and the increase of the driving force of 
formation of locally dense packed states, which progressively favor nucleation of the "liquid" 
(colloid-rich) phase|4|. As a result, the establishment of conditions such that a connected 
percolating structure able to sustain stress survives for times longer than the observation 
time is often preempted by phase separation. For spherical attractive interaction potentials, 
it has been found that phase separation is always dominant at low densities and arrested 
states at low packing fractions (f) can be reached only in non- homogeneous state, whose 
morphology is determined by the phase separation process 13|, [ij, |l5|, |l6|, |l7|, |l8|, |l9| . 



A possible mechanism to generate low packing fraction arrested states in the absence of 



2l|. 



phase separation has been recently proposed and supported by off-lattice simulations 
It has been shown that the region in the T — cf) plane in which a two-phase coexistence is 
thermodynamically preferred as compared to the homogeneous fluid state can be progres- 
sively reduced (both in T and 0) by decreasing the maximum number of possible pair-wise 
bonded interactions. In an equivalent language, the reduction — at fixed interaction poten- 
tial range — of the colloid surface allowing for inter-particle bonding progressively destabilize 
phase separation^, According to these ideas, colloidal particles with a limited number 
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of attractive spots are the best candidates for gel formation. Interestingly enough, these 
ideas carry also to the case of protein solutions, where the character of the amino acids on 
the surface of the protein controls the strength and the directionality of the inter-protein 
interaction P 

In this article, we present a lattice gas model, which we solve in the Bragg- Williams 
mean field approximation and in the Bethe-Peierls quasi- chemical approximation and which 
allow us to visualize in a clear way the relation between the limitation of the maximum 
number of interaction and the amplitude of the phase separated region in the T-0 phase 
diagram. We complement these calculations with Monte Carlo simulations of the liquid-gas 
phase coexistence, to assess the reliability of the mean field solutions. The reported results 
confirm that the reduction of the maximum number of interactions is indeed an efficient 
mechanism for generating thermodynamically stable states at extremely low temperatures. 

II. THE MODEL 

The system we consider is a nearest neighbor lattice gas, with occupancy variables Ui at 
each node of a lattice with 7 nearest neighbors. 

The Hamiltonian for the system is written, in the usual form, as 

<ij> 

but with the * indicating that the sum is only over such bond configurations that have 
at each vertex a maximum of 7^ bonds. 

For the purposes of mean field theory, we can write this as 

i j 

where the first sum i runs over all lattice sites, and the second sum j runs over all 7 
nearest neighbors of site i, and with 

f{x) = X, X <-fm (3) 
= -fm, X> 7m (4) 

where 7^ denotes the maximum number of interactions, or valency, allowed for a particle 
at any site. When 7^ = 7, one recovers the simple lattice gas. In the following, we consider 



the behavior of the system when the parameter 7^ is varied between the hmits 7m = 
(when one has a non-interacting paramagnetic lattice gas) and 7to = 7 (when one has the 
simple nearest neighbor lattice gas). 



III. MEAN FIELD APPROXIMATION 



In order to perform the mean field calculation, we approximate the function / as 

/(a:^) = 7mtanh(x/7m) (5) 

which has the desired linear behavior for small x and a constant value of 7^ at large x. 
Further, the mean field approximation amounts to writing Uj = < n > where < n > 
is the average occupancy, given by < n >= rii/N where N is the total number of sites. 
With this approximation, the thermodynamic potential Q in the grand canonical ensemble 
may be written as 



e 7 < 77- > 
Q/N = --7m <n> tanh(-^ )-// <n> ^ksT [< n > log < n > +(1- <n>) log(l- < n >)] 

2 7m 

(6) 

where 11 is the chemical potential, T is the temperature, ks is the Boltzmann's constant 
and < n > is the occupancy, given n and T, that minimizes Vt. Imposing the minimization 
condition. 



d <n> 

allows one to eliminate /j, and write the thermodynamic potential Q in terms of T and 
< n >. Further, since Q = —PV = —PN, (where P is the pressure, and the system volume 
for the lattice gas V — N the number of sites) one obtains in this way the equation of state 
of the system: 

P = < n sech^^^^^) - keTlogil- < n >) (7) 

7™ 

The condition for the critical point is given by 

o < n > 
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FIG. 1: Mean field critical density (top panel) and critical temperature (lower panel) as a function 
of parameter 7m/7- Calculated for 7 = 12 and e = 1. 

and 



0. 



(9 < n >^ 

From these, one obtains the critical density and the critical temperature, which depend 
on the parameter 7^/7 and are plotted in Fig. 1. From the figure, it is clear that both the 
critical density and the critical temperature go to zero as the variable 7^/7 goes to zero. 

We next calculate the spinodals that demarcate the region of instabihty in the (< n > 
,T) and in the (< n >,P) phase diagram which are shown in Fig. 2. Consistently with 
expectations based on the behavior of the critical density and temperature, the region of 
instability also shrinks with the reduction in the valency. 
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<n> 



FIG. 2: The mean-field spinodals in the (< n >, T) phase diagram (top panel) and in the (< n >,P) 
phase diagram (lower panel) for different values of the parameter 7m- Calculated for 7 = 12 and 



IV. BETHE-PEIERLS APPROXIMATION 

As an improvement over the mean field approximation, we next consider the Bethe- 
Peierls approximation, which in the case of the sirnple lattice gas is equivalent to the Bethe 
lattice. As in the standard treatment (see, e. g. |2J]), we consider a "center" site, and its 
surrounding neighbors. We write the probability that the center site is occupied, with n 
nearest neighbors also being occupied, as 



e = 1. 



P{l,n) = i exp(/?/i)C^ exp(/3en) exp(/5/in)z" ^ < 7, 
= iQ exp(/3e7m) exp(/5/in)z" n > 7, 



m 



m 



(8) 
(9) 
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where (5 = I/ZcbT, q is the normahzation, and z is introduced to account for the influence 
of the rest of the lattice. The rest of the terms explicitly account for the interactions and 
the chemical potential corresponding to the number of occupied sites present. Similarly, the 
probability that the center site is not occupied, with n neighbors occupied, is 

P(0, n) = exp(/3/xn)^" (10) 
The normalization factor q is given by the condition 

7 

J][P(l,n) + P(0,n)] = l (11) 

n=0 

The average occupancy number < no > for the the center site is the probability that the 
center is occupied regardless of the occupancy of the surrounding sites, that is: 

7 

<no>=^P(l,n). (12) 

n=0 

The average occupancy number < Ug > for a given neighbor site is given by 

1 ^ 

< rie >= -V n[P(l,n) +P(0,n)]. (13) 

' n=0 

The unknown parameter z is now determined by the condition that these two occupancy 
number, < tiq > and < ng > are both equal to < n >. Defining a function 

g{z, fj,, T) =< no > - < Ue >, 

the solution is found by the condition g{z, n, T) = 0. Determining z in this way, and using 
the expression for < no > above, one obtains the probability of occupation as a funciton of 
T and /i. In addition, imposing the condition that the first and second derivatives of g also 
vanish yeilds the condition for the critical point, since this vanishing of the first and second 
derivatives marks the change from the high temperature regime where only one solution 
exists, to the presence of more than one solution at lower temperatures. In addition, the 
condition that the first derivative vanishes is used to locate the spinodal points. Figure 3 
shows the dependence of the critical density, temperature, and the chemical potential on the 
parameter which is varied in integer steps for 7 = 12. The case of 7^ = 7 reproduces 
the standard result for this approximation. It is seen that the critical temperature decreases 
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FIG. 3: Critical density and critical temperature as a function of parameter 7^ for the Bethe-Peierls 
approximation . 

very slowly at first as 7^ is reduced, but goes to zero for 7^ = 1. Interestingly, the case 
of 7m = 2 displays a finite critical temperature, which is somewhat surprising, since in this 
case the system can form only two bonds at each site, and one expects this to be equivalent 
to a one dimensional system with a vanishing critical temperature. 

We next calculate the spinodals that demarcate the region of instability in the {< n >,T) 
phase diagram for 7m = 3, which is shown in Fig. 4. As in the mean field case the region of 
instability also shrinks with the reduction in the valency. 

V. SIMULATIONS 

In order to estimate the accuracy of the approximations above, we evaluate the critical 
parameters and the co-existence lines for different values of 7m using computer simulations, 
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FIG. 4: The Bethe-Peierls spinodal line evaluated for the specific case of 7 = 12 and 7m = 3 . 



emplying the histogram reweig hting technique H H Q • 

In this method, histograms of 

sampled energies and number densities from simulations at different external parameters are 
used to estimate composite probability distributions, which in turn may be used to obtain 
information on phase equilibrium. To this end, we perform Monte Carlo simulations in the 
constant temperature and constant chemical potential ensemble for a lattice of 256 sites 
arranged in an FCC lattice. The energy of the system depends not only on which sites are 
occupied, but also which bonds exist in the system. Therefore, care must be taken that the 
formation and deletion of bonds obeys detailed balance. In the grand canonical simulations, 
it is sufficient, when a particle is inserted, to pick any of the possible bond configurations at 
random, with equal probability. At each simulated temperature and chemical potential, runs 
of length 2.5 million Monte Carlo cycles (MCS). The energies and number of particles every 
100 MCS, after an equilibration of 0.5 million MCS, are used to build histograms fi{N,E) 
(where N is the number of particles and E is the energ y), where i = 1 . . . R indexes the 
different runs at different /x and T values. As described in j^, , the composite probability 



distribtion r(A^, E, fj,, (5) is obtained with 



r(iv,E,/x,/3) 



(14) 



where is the total number of observations for each run, and the constants Cj are obtained 
iteratively from the relationship 



Once the above equations are run to convergence, the probability distribution in particle 
number can be obtained by summing F over the energies. Phase co-existence at any given 
temperature is obtained by requiring that the distribution with respect to density has two 
distinct peaks, of equal height. Figure 5 shows the coexistence curves obtained for the 
cases jrn — 12, 9, 6. Studying phase coexistence at lower values of is severely hampered 
by the extremely slow equilibration in the temperature range of interest and has not been 
attempted. In each of the cases studied, evaluating the coexisting phase densities is limited 
by the ability resolve the coexisting density peaks, which is difficult as the critical point 
is approached. Nevertheless, it is clear from the data shown that the simulations confirm 
the trend seen in the calculations above, that the critical tempeature, density and (minus) 
the chemical potential decrease as 7^ decreases. Figure 6 shows the critical parameters 
in comparison with the Bethe-Peierls calculations. It is seen that the critical paramters 
drop more rapidly than suggested by the Bethe-Peierls calculations. Unlike the 7^ = 12 
case, where the coexistence chemical potential and the mean of the coexisting densities are 
constant below the critical temperature, both these quantities increase as one moves to lower 
temperatures below the critical temperature. 

VI. CONCLUSIONS 

We have presented results for lattice gas models with the imposition of a constraint on the 
maximum number of bonds that particles can participate in. Mean field (Bragg- Williams 
and Bethe-Peierls) calculations and computer simulations using the histogram reweighting 
technique, show that the critical temperature and density decrease as the maximum number 
of bonds allowed for a given particle is reduced. Also, we find that the density range of the 
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FIG. 5: The coexistence curves for 7^ = 12, 9, 6 from simulations, showing that the coexistence 
region shrinks as 7^ decreases. 

coexistence region, where the hquid and gas phases coexist, shrinks as the maximum valency 
of the particles is reduced. These results are consistent with the results that have been 
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2l\ . It is thus confirmed that 



obtained earlier for similarly defined continuum models 
valency reduction is effective in opening a large region of densities in which low temperature 
disordered states can be accessed in equilibrium. Studies of models with controlled valency, 
both on the lattice and off-lattice, may make it possible to disentangle gelation from phase 
separation. 
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FIG. 6: Comparison of the critical parameters from simulation with those of the Bethe-Peierls 
approximation, indicating that the critical parameters decrease more rapidly than suggested by 
the Bethe-Peierls approximation. 
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